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NUMERICAL  SOLUTIONS  OF  POLYNOMIAL  EQUATIONS 

ABSTRACT 

This  report  describes  the  FORTRAN  Subroutine  POLYR  and  a  related 
complete  program  BRL-RSSR  for  finding  all  roots  (real  and  complex)  of  a 
real  polynomial  equation 

N 

P(x)  -  V  A.XK_1  =  0. 

i=0 

v 

The  method  which  is  used  combines  the  root  squaring  and  the  subresultant 
(extracting  quadratic  factors)  procedures  for  calculating  all  roots,  e\  n 
multiple  roots,  of  real  polynomials.  Reconstruction  of  the  coefficients 
from  the  product 

N 

n  (x  -  Root. ) , 

1=1 

which  is  included,  serves  as  a  means  of  checking. 

The  described  subroutine  and  program  were  adapted  from  a  similar 
program  in  double  precision  described  in  Report  ALL  6967,  Argonne 
National  Laboratory  by  Erwin  H.  Bareiss  and  Ronald  Hamelink. 
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ROOTS  OF  POLYNOMIALS 


Subro  tine  POLYR  and  BRL-RSSR  Program 

General  Description.  A  FORTRAN  subroutine  and  a  related 
complete  program  for  finding  all  roots  {real  and  complex)  of  a  real 
polynomial  equation 

N 

P(x)  =  )  A^vN  1  _  u 
i=0 

are  available  for  use  on  the  BRLESC  computer. 

Method.  The  rest  squaring  and  subresultant  (extracting  quadratic 
factors)  procedure  "RSSR"  for  finding  real,  complex  and  repeated 
roots  of  real  polynomials  and  checking  them  by  reconstructing  the 
coefficients  from  the  product 

N 

n  (x-ROOT.) 
i=l 

is  described  in  detail  in 

Report  ANL  6987 

Argonne  National  Laboratory 

by  Edwin  H.  Bareiss  and  Ronald  Hamelink  . 

Authors.  The  routine  under  the  name  RSSR  was  written  in 
FORTRAN  by  the  above  authors  for  the  IBM-704  and  CDC-3600 
computers  using  double  precision  arithmetic.  The  modifications 
and  additions  necessary  to  make  it  run  on  BRLESC  in  single 
precision  were  performed  by  Henry  Wisniewski  of  the  Computing 
Laboratory  of  BRL  . 

Modifications  to  the  Subroutine  POLYR.  When  compared  with 
the  original  FORTR  I  ANL  subroutine  (pp  49-73)  the  BRL  subroutine 
POLYR  differs  as  follows: 
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a)  The  Main  Program,  which  does  input  and  output  and  calls 
RSSR,  has  been  entirely  replaced  (see  tabulation). 

b)  Subroutine  RSSR 
No  changes. 

c)  Subroutine  ROOTSQ 
Changes: 

LINE  STATEMENT 

1+6  IF  (KM)  11,  4,  11 

11  DO  4  L=l,  K  M 

d)  Subroutine  REALROOT 
rv.anges: 

LINE  STATEMENT 

11  +  1  T=2**M*I 

Reason:  Mixed  type  arithmetic  is  not  allowed  in  BRLESC 
FORTRAN. 

e)  Subroutine  COMPROOT 
No  changes. 

f)  Subroutine  TEST 
Changes : 


LINE 

STATEMENT 

11 

IF  (IB(3)  -  74) 

100, 

100, 

12 

100 

IF  (I B( 3 )  +  74) 

12, 

101, 

101 

18 

IF  (IB(2)  -  74) 

102, 

102, 

12 

103 

IF  (IB(2)  +  74) 

12, 

103, 

103 

103+1 

IF  (IB(3)  -  74) 

104, 

104, 

12 

104 

IF  (IB(3)  +  74) 

12, 

105, 

105 

g,h,i)  Subroutines  SUBRES,  RECON,  QUADRIV 
No  changes. 

j)  Subroutine  DOUBLOG 
Changes: 

LINE  STATEMENT 

1  PRINT  2 

3  To  =  lx 

3+1  y  =  4LOG(x)  +  To*ALOG(T) 


b 


Rea  son;  Mixed  type  arithmetic  is  not  allowed. 


k.)  Subroutine  DOUBLEXP 
Cnanges : 

LINE  STATEMENT 

First  Z  =  EXP(X)  $  IZ  =  0 

1)  Subroutine  ADD 
Changes : 

LINE  STATEMENT 

10  IF  (IDIFF)  11,  12,  11 

12  DO  111=1,  IDIFF 

m,n)  Subroutines  SUBTRACT  and  SCALE 
No  changes. 

o)  Subroutine  UNSCALE 
Changes : 

LINE  STATEMENT 

First  +  1  IF (IX  +  84)  3,  4,  4 

4  IF(IX  -  84)  5,  =>,  0 

Reason:  Max;mum  size  of  IX  on  BRLLSC  is  84 
o  X  --  1. 0E  f  153 

Reason:  Maximum  number  on  I3RLESC  is  1.  34078EIS4 
0*2  PRINT  7 

7  FORMAT  (25H0EXP.  OVERFLOW  [N  UNSCALE) 

Operating  Instructions  for  the  Subroutine  PGLYR,  They  are 
replaced  by: 

Entry  is  made  by  the  statement 

CALL  POLYR  (N.  COEFF.  ROOTS.  D) 

v.  he  re 

N  =  Degree  of  the  polynomial.  (A  lixed  point  integer.  ) 
COEFF  =  One-dimens lonal  array  of  length  (N+l)  which  contains 
the  locations  of  the  coefficients. 

ROOTS  =  Two-dimensional  array  of  length  (2xN)  which  contains 
the  roots. 


D  -  One-dimensional  array  of  length  (N+l)  which  contains 
the  reconstructed  polynomial  from  the  roots  and  the 
first  coefficients. 

The  coefficient  of  the  h.^  iest  degree  is  stored  first, 
followed  by  the  remaining  coefficients  in  decreasing  order  of  degree. 

Storage 

COEFF  (1)  A 

o 

COEFF  (2)  Aj 

COEFF  (3)  A2 

COEFF  (N+l)  A.. 


The  results  are: 

ROOTS  (1,1)  Real  Part 

ROOTS  (1,2)  Real  Part 

ROOTS  (1.3)  Real  Part 


ROOTS  (l.N)  Real  Part 


ROOTS  (2,  1)  Imaginary  Part 

ROOTS  (2,  2)  Imaginary  Part 

ROOTS  <2..  3)  Imaginary  Part 


ROOTS  (2,  N)  Imaginary  Part 


The  coefficients  of  the  reconstructed  polynomial. 
D(l)  An 


D(N+  1 )  AN 

This  subroutine  is  available  on  cards  only.  (See  Appendix  A) 
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Modifications  to  the  Program  RSSR.  The  BRL-RSSR  Program 
contains  the  same  modifications  as  POLYR  with  the  exception  of  a 
Main  Program  which  is  retained  as  shown  in  Report  ANL  b9B7 
with  the  following  changes: 


LINE  STATEMENT 

2-rl  and  3  Cancelled 


4  PRINT  2,  TITLE 

4+3  Cancelled 

10+1 


14+1  through  18 

19  +  2 


10+  \  END 

31  through  37*2  Cancelled 


Operating  Instructions  for  the  BRL-RSSR  Piogram.  As  in 
Report  ANL  o987,p.  73,  with  the  following  changes: 

a.  No  change 

b.  No  change 

c.  Cancelled 

d.  The  first  line  is  changed  as  to  Hows: 

A  .  Format  .  £  18.  10  • 

The  remaining  lines  are  not  changed.  <See  Appendix  B) 

Example 

The  input  and  output  ot  the  BRLESC  computations  of  the  roots  oi  the 
very  badly  conditioned  polynomial 

^  4  3 

P(x)  z  7 0 *>(.' .  x  -  MOT 2.x  *  2  3t>x  ♦  <>t?  i  if> 0,  64x 

> 

-  2  i  0  1 0 H 4 .  L x "  4  x  -  104KMv>.  -  0 


are  given  in  Appendix  C. 

A',  ^liability.  Copies  of  the  punched  iard  decks  .or  POLVR  and 
RSSR  iii.iv  be  obtained  from  Systems  Progr  .milling  Branch.  Computer 
S  ;p!>ort  D.  vision. 
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APPENDIX  A 


SUBROUTINE  POLYRIN,COEFr,ROGlS,i))  l 

DIMENSION  A(  SI  ,31  ,|  A(  si  ,  J)  ,  ROD  f  S(  1  ,N  )  ,  I)  <  1 )  ,  COEFF  <  1 )  2 

TYPE  1  N  T  r  G  f  K  UEORFF  3 

DESK PH  =  N  4 

N1*DEGREF+1  5 

M  =  iu  6 

MMAX=15  7 

DELTA=0.C001  8 

EPS ILCN-j.OvOGOl  9 

DC  14  I  - 1 »N1  10 

All  ,l)=COFFFU)  11 

IMI.D--0  12 

CALL  SCALE IA( I ,11 ,1 A( I ,1 ) I  13 

14  CONTINUE  14 

CALL  RSSKIA.I A, R DOTS , DEGREE . M , MK AX , DEL T A , EPSILON, 0)  15 

1 F ( N 1- I  DEGREE ♦! ) 129,29,26  16 

26  RETURN  17 

29  PRINT  3o  18 

30  ECRHATI 21H0SGME  ROOTS  NOT  FOUND)  19 

RETURN  20 

END  21 


SUBROUTINE  R  S  SR  ( A  ,  I  A , R 00 TS , DE GRE E , M , MM A  X , DEL T A, E PS  I  LON , 0 ) 
DIMENSION  A! 51,3) , I A(  51 ,3) , ROOTS! 2, SC ) ,D(51 ) .RCMODl 50 ) , MROMOU ( 50  )  , 
1  NON R  T I  5  J ) , MNQNRT ( 50 ) 

TYPE  INTEGER  DEGREE 
N- DEGREE 
IF  (N)  1,1,3 

1  dec.ree  =  ncur 

2  RETURN 

3  N 1  *  N  *  1 
N2  =  NU1 

00  5  1  =1  , N 
K=N2-I 

IF  ( AIK.l  )  )  6,4,6 

4  J  *  N 1  *  I 

ROOTS! 1 ,J) =0.0 
R0CTSI2, J)*C.  0 

5  CONTINUE 
DEGRCE=C 
GO  TO  2 

6  N1*K 
N  =  K-l 
NCUR=N 
NL=N 

7  CALL  ROOTSO  I  A , l A , NCUR, M I 

CALL  REALROOT  (A, IA, M.NCUR, DELTA, EPSILON, ROMOD,MROMOD, 

1  NONA T.MNCNRT.NCC, ROOTS) 

IF  I  NO  0 )  12,12,4 
6  N1«NCUR+1 

CALL  COMPROCT  (  A  ,  |  A  ,  •<  CMOD  ,ROO  T  S  ,  M  ,  MNONRT  ,NDNRT  , 

1MKCM0D.NC0, DELTA, EPSILON, NCUR) 

IF  (NCUR)  12, .2, 9 
9  IF  (NL-NCUR)  11,11,10 

10  Nl  *N'„UR 
GO  TC  7 

11  M  *  M ♦ 1 

IF  IMMAX-M)  1,7,7 

12  CALL  RECON  (ROOTS, All, l) ,IA(l,li,D. DEGREE) 

GO  TO  1 


U 


END 

SURROUTINF  ROOTSO  I  A , l A , NCUR, MM  I 
DIMENSION  A {  5i ,3 ) , I A( 51 »  3  ) 

N1  =NCUR ♦ I 
DO  1  J-  1 » N1 
A(J»2)=A(J»i) 

IA(J»2)-IA(Jfl) 

A I  J ,  3  )  *  G  •  0 
I A  (  J, 3 ) =C 

1  CGNTINUC 

OG  9  M*  1 , MM 
DO  6  J=l,Nl 
Kl *Nl - J 
K2  =  J- 1 

KM=XMIN0F(K4,K2) 

IE  (  KM  J  11,  •♦,11 
11  00  4  1=1, KM 
LR*XMC0F(L,2) 

JL* J-L 

JLP-J+L 

IF  (LR)  2,2,3 

2  X*A(vIL,2)*A(JLP»2) 

I <  = I  ft ( JL,2  » ♦ I  A  < JLP,2) 
call  SCALF  ( X ,  I X ) 

CALL  ADO  (A( J,3) ,IA( J,3) ,X, IX, A( J»3) ♦ l A( J,3)  ) 

GO  TO  4 

3  X»A< JL,c)*A( JLP.2) 

IX=IA(JL»2)*IA(JLP»2) 

CALL  SCALE  (X,IX) 

CALL  SUBTRACT  ( A r J , 3 ) , I  A ( J , 3 ) , X , I  X , A ( J , 3 ) , I  A ( J , 3  I  ) 

4  CONTINUE 

A ( J , 3 l =2 • 0*A  ( J  ,  ?  ) 

CALL  SCAL^  (A(J,3),IACJ,3n 
X  =  A( J,2  1**2 
IX=IA(J,2)*IA(J,2) 

CALL  SCALC  (X,IX) 

CALL  ADD  ( A( J,3) , I A< J , 3  » , X, IX, A< J,3> , I  A  C  J , 3)  ) 

JR=XMCDF  ( J , 2 ) 

IF  (JR)  5,5,6 

5  A(J,3)*-A(J,3) 

6  CONTINUE 

IF  (MM-M)  9,9,7 

7  DC  «  J-l.Nl 
A(J,2)=ACJ,3) 

IA(J,2)=IA(J,3) 

A(J,3)*0»0 

I A I J , 3 ) =0 

8  CONTINUE 

9  CONTINUE 
RETURN 
END 

SUBROUTINE  REALROOT  U , l A ,M ,NCUR ,OELT A, EPS  I  LON, ROMQD, MRGMQU, 
lNONRT,MNONRT,NCO, ROOTS) 

DIMENSION  Ai51*3)»IA(51»3), R00TSI2, 50 >» ROMO 0(50), MROMODI 50 ) , 
1N0NRTI 50) ,MNONRT (50) , RAT  10(51 1 , I P I V I  5 1 ) ♦ AREOI 50 ) , I AREDI 50 ) 
RATICI1)=1.0 
DC  4  1-2, NCUR 
I  i-XMCDF (1,2) 


} 
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IF  ( A  C I , 3  )  >  2,1,2 

1  RATIO! t)=C.O 
GC  TO  4 

2  T  =  A ( I,2)*A( 1,2) 

I  T*  I A  t  I  ,?>*IA!  I,  2) 

CALL  SC4LL  IT, IT) 

T=T/A< 1,3) 

I T= I T- I  A  I  1,3) 

IF  (IT-2)  50,50,1 

50  IF  <IT*2)  1,51,51 

51  CALL  UNSCALE  I T, I T) 

RATIO! I  )=T 

IF  (ID  3,3,4 

3  RATIO! 1 ) =-RATI 0 (  !  ) 

4  CONTINUE 
RATI0(NCUR+1)=1,0 
IP  I V ( 1 )  =  I 
IPIV(NCUR*1 )=1 

00  7  l =  2  » NC  UR 
K=ABSF(R4TIC( I ) - 1 . 0 ) 

IF  l X-DFLTA )  5,6,5 

5  IP  I  V! I ) =1 
GC  TO  7 

6  I P  I  V  (  I  )  =  C 

7  CONTINUE 
INCUR  1 =NCUR ♦  1 

II  =  0 
MUL  T  =  0 
1  =  1 

14  =  1 

8  11*11*1 
12=11*1 

MUL  T=MULT*  1 
IF  ( IPI V! 12)  )  8,8,9 
9  RCM0D(I4)=a!I2,3)/A(! ,3) 

I KOMOD= I  A ( I  2  »  3 ) -  I A (  1,3) 

CALL  SCALE  !  ROMOO !  1 4 )  ,  I  RQMOU) 

IF  ! RCMOD ( I  4 ) )  10,11,11 

10  ROMOO!  I  4  )  =-ROMOL) !  14) 

11  CALL  OOURLOG  ( R 0 MQD ( I  4 ) , I ROMOO , XN , I XN ) 

T  =  2  **M*  1 1 

XN=XN/T 

CALL  SCALE  (XN.IXN) 

CALL  DOURLEXP  !  XN ,  E XN , ROMOD ! I  4 ) , I ROMOD ) 
I F ( I RCMOO- 74  )  100,100,101 

100  IF! IROM0D+74  )  101,101,102 

101  ROMOD! I4)=0.0 
IROMCD*u 

GO  TO  103 

102  CALL  UNSCALE  ( ROMOO ! I  4 ) , I ROMOO  ) 

103  MROMOD ! !4)=MULT 

IF  ! NCUR* 1-  I  2  )  13,13,12 

12  1=12 
14=14*1 
MUL  T  =  0 
11*0 

GO  TO  8 
13  y=c.c 
NC0=O 

00  22  1=1,14 
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KL=  I  <»♦  1- I 
W  =  *ROMOD ( KL  } 

[ 5*MRCM0D ( KL ) 

00  20  J  =  1 »  I  5 
J*J 

14  CALL  TEST  ( A,  I A , W. Q,NCUR , RQMOD ( KL ) » EPS  I  LON, K \ 

IF  (K)  17,17,15 

15  ROOTS! 1 ,MCUR : 

R00TS<2  »NCUR)=Q.O 
ARELH  1 )=A<l,i  ) 

I AR6D ( 1  I  =  I  A ( 1  , 1 ) 

DO  16  L=2,NCUR 
Y=ARE0(L-1)*W 
I Y» I AREO ( l -1 ) 

CALL  SCALE  <Y,[Y) 

CALL  SUBTRACT  I A ( L , U , I  A ( L , 1 ) , Y , I Y , ARED ( L ) , I AREOI L I) 

A( L , 1 ) sARCO (L  ) 

I A ( L , 1 )  =  I AREOI L ) 

16  CONTINUE 
GO  TO  19 

17  IF  (Ml  18,21,21 

18  W=-W 

GO  TO  14 

19  NCUR=NCUR-l 

20  CONTINUE 
GO  TO  22 

21  NCONCO+1 
NONRT ( NCC ) =KL 
MN0NKT(NCG)=15*1-J 

22  CONTINUE 
RETURN 
END 

SUBROUTINE  COMPROOT  ( A, I  A , ROMOD , ROC TS , M, MNON "T , NONRT , 

1 MR CMOD,NCO, DELTA, EPSILON, NC UR  I 

DIMENSION  A(51,3),lA(5l,3),R0M0O(5O),R00TS(2*5O),SR'e'l,3),ISR(51,3 
1 1.SR0M0UI50)  ,SROOTSI2,5Q)  ,  MNONRT I  50  »  ,  NONRT(  5O),MSR0MOD(  50, 
2NSCNRTI49) , MSNOR  T  ( 49 ) ,MKOMOO(50) ,D(2),R(2),RC49) 

00  23  1=1, NCO 
JA=NONR  T ( I ) 

U*MNONRT(  I  ) 

11*11/2 
IF  (ID  1,1,2 

1  11*1 

2  IF  I ROMODI JAI >  3,23,3 

3  Q*R0MCO(JA) 

DO  22  J» 1 , 1 1 

300  CALL  SUhRFS  ( A , I  A , NCUR, SR , I  SR , 0) 

IF  (NCUR-4)  4,100,5 

100  NSCUR *2 
GO  TO  101 

4  NSCUR*1 
J2*l 

GO  TC  6 

5  NSCUR *NCUR- 3 

101  J2*NSCUR 

6  LL  *N SC UR* 1 

IF  ( NSCUR- 1 )  7,7,9 

7  IF  I  SR  1 1 , 1)  )  8,10,8 

8  X*  SK I  1 , D  «  I  X*f  SR  1 1, 1 ) 
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V  =  SPM  £  ,  1  >  V  I  V=  I  S=i  I  2  *  1 » 
call  UNSCALF  IX, t  X  > 

CALL  C  ;  SC  ALE  ( Y, \ V) 

SRC  0  T  S ( 1,1 ) =  -Y/X 
NSCUR=0 

go  rc  n 

9  CALL  ROJTSO  ( SR , I  SR ,NSCUR  ,M > 

CALL  REALROCT  ( SRt l SR , M , NSC UR , UtL T A , EPS  I  LON, SRCMOD, MSROMCD, 
1NS0NRT,MSN0RT,MSCD, SROOTS) 

IE  ( J2-NSCUK )  1C,  1C,  11 

10  SRCOTSI 1, J2  1  =C  . 0 

11  SROGTSI 1 ,J2 )= SROOTS tl ,J2 )*RGMOO( JA) 

T=ROMODUAI+RCMOD(  JA) 

IE  ( SROUl S( 1  *  J2 ) -T )  12,16,16 

12  W-SROCTSI 1, J2) 

WE  =  R0MGD(  JA  )*ROMC'D(  JA  ) 

CALL  TEST  ( A, I  A, W,WC ,NCUR,ROMOO( JA ) , EPS IL0N»K  I 
IF  (K)  15, IS, 13 

13  ROOTS! 1 ,NCUR) =-W/’.0 
T=4.0*WE 

U=W*W 
T  =  T-U 

IF  (T)  201,201,202 

201  T=-T 

U  =  l)SQRT  (  r ) 

ROOTS! 1 »NCUR) =ROCTS ( 1 *NCUR)-U/2.0 
ROOTS! l,NCUR-i )=-U-U 1/2.0 
ROC  TS ( 2  » NCUR ) =0, 0 
ROOTS! 2»NCUR-i 1=0.0 
GO  TO 

202  U=DSQRT ! T ) 

203  ROOTS! 2,MCURI=U/2.C 

ROOTS! 1,NCUR-1 I =R00 TS ! I , NCUR ) 

ROC  TS! 2,NCUR-i ) =-ROCTS( 2 ,NCUR) 

204  D! 1 )=W 
0(2  >=WE 

CALL  OUADIV  ( NCUR, A , I  A, R , 0, B I 
JX=NCUR-1 
00  14  JY=1 , JX 
A!  J Y , 1 )=R{ JY) 

I  A ( J Y , 1 ) =C 

CALL  SCALE  ( A ( JY , 1 ) , I  A ( J Y  ,  1  )> 

14  CGNTINUE 
NCHR=NCUR-2 
GO  TO  22 

15  W=-W 

CALL  TEST  ( A , I  A , W, WE , NCUR , RuMOU! J A ) , EPS  I  LON, K I 
IF  (K)  16,16,13 

16  IF  U2-INSCUR  +  1))  23,17,19 

17  IF  ( J  2  —  1  >  23,23,18 

18  J2= J2- 1 

SROOTS 1 1 , J2 I *0.0 
GO  TC  12 

19  IF  (SR0uTS(l,J2)-SR00TS(l,J2-m  20,21,20 

20  J2* J2-1 
GO  TO  11 

21  J  2  =  J  2 - 1 
GO  TO  16 

22  CONTINUE 

23  CONTINUE 


% 
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RETURN 

END 


SUBROUTINE  TEST  ( A , I  A , * , 0 , N , ROMGD , EPS  I LCN , K ) 

DIMENSION  A(51,3),IA(51,3),B(3),IB(3>,T(2l,E(2),C(51i 
B ( 1  )=0.0$l  X  =  0UW  =  0 
I B ( 1 1 =0 
R (  2 ) =A ( 1 , 1 ) 

IB(2)=IA(1,1) 

DO  2  1  =  1,  N 
X  =  W*R ( 2  ) 

I X*  I R ( 2  ) 

CALL  SCALE  IX,IXI 
Y=Q*B ( 1 ) 

I Y  *  I  fl  ( 1  ) 

CALL  SCALE  l Y , I Y ) 

CALL  ADO  (X,IX,Y,IY,Z,IZ) 

CALL  SUBTRACT  ( A (I U , 1 ) , I A ( I  +  1 , 1 ) , Z , I Z , B ( 3 ) , IB ( 3 )) 

IF  (N-I)  2,2,1 

1  B  ( 1 ) =B ( 2  ) 

IB  I  2  )  =  I  0  (  2  I 
B I  2  I =P ( 3  I 
IHI 2  5  =  I B ( 3 ) 

2  CONTINUE 

3  K9UNT*1 
CEPSIL*EPSIL0N 
T(l)*C.OiT(2)=0.0 
N1=N>1 

X=2 • 0*EPSl LON 

Y=X*RCMGD 

6(11 =ROMOD+ Y 

E(2)=ROMOO*CEPSIL*ROMGD 

DO  7  1*1, Nl 

IF  (All,!))  4,5,5 

4  C< l !=-A< I ,i)*IC=l A( I ,1) 

GO  TO  6 

5  C( I  )  = A (  I,1)$IC  =  IA(I,1) 

6  CALL  UNSCALE  (C( l  l,IC) 

T(ll*T(i)*Flil*C(ll 
T(2)*T(Z)*E(2)*C(  I  ) 

7  CONTINUE 
DIF*T(1)-T(2) 

8  IF  (01  18,9, IH 

9  IF  <  B( 3  1  )  10,11,11 

10  B( 3  )  *-B ( 3  I 

11  IF  ( I B ( 3 ) *74  )  100,100,12 

100  IF  ( IB ( 31*74  )  12,101,101 

101  CALL  UNSCALE  (8(31,18(3)1 
IF  ( 0  I F-B ( 3 )  )  12,12,17 

12  K«0 

IE  (KOUNT-2)  13,16,16 

13  IF  (Q)  15,14,15 

14  IF  (W)  15,16,i6 

15  SENSE  LIGHT  2 
KGUN  T*KOUNT ♦ 1 

16  RETURN 

17  K*  1 

GO  TO  16 

18  IF  (18(2 )-74  )  102,102,12 

102  IF  ( 1 8 ( 2 ) ♦ 74  )  12,103,103 
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103  CALL  UN  SC ALL  ( 8 ( 2 )  ,  I R C 2  )  ) 

IF  (  IP(  3)-?4  )  104,104,12 

104  IF  ( l A ( 3 )  +  74  )  12,105,105 

105  CALL  UN  SC  ALL  (HO), 10(3)) 

X=Q*R ( 2  )*R( I  ) 

Y=W*R(2)*R<3) 

Z=B ( 3 ) *B ( 3 ) 

Y=  X~  Y  +  Z 

IF  (Y)  10,17,20 

19  Y  =  -Y 

20  0 1 F  =  0 l F * 0  I  F 

IF  (0IF-Y)  12,17,17 
END 

SURK  OUT i NF  SURRFS  (A, I  A , N , SR , I  SR, ROMO 0) 

DIMENSION  A C51,3),IA( 51,3 ) , SR ( 51 , 3 ) , I  SR ( 51,3), CC51),B( 50,3) 

N1=N*1 

T*1.C 

DO  1  I =1  ,N 
J=N 1- F 
T  =  T  *R  CMOD 

C( J)=A( J,1)*T 
1C-IA(J,  1) 

CALL  UNSCALE  (C(J»,IC) 

1  CONTINUE 

C  (  N  l  )  =  A  (  N1  ,  1  ) 

IC=IA(Nl,i> 

CALL  UNSCALE  (C(Ni)»IC) 

IF  (N-2)  12,12,2 

2  N2  =  N-2 

DO  3  I = i , N2 
a(  i , l )  =  c . o 
R(  I  ,  2  )  =  0 . 0 

3  CONTINUE 
1=2 

B(1,2)=C(1> 

4  8 ( 1 , 3 ) =C ( I >-~ (1,1) 

DO  5  J=2,N2 

B( J,3)=-G( J-l,2)-R( J,l) 

5  CONTINUE 

IF  ( N- ( 3 ♦ l I i  100,6,6 

6  I  =  I  ♦  1 

DO  7  J  =  1 ,  N 2 
B(J,1)*R(J,2) 

B(J,2)=H(J,3) 

7  CONTINUt 
GC  TO  4 

100  IF  (N-4)  14,101,0 

101  IF  (N-(2*U)  104, 102,102 

102  1=1*1 

00  103  J  =  1 » 2 
B(J»l)*fl(J,2) 
fll J,2)=R( J,3) 

103  CONTINUE 
GC  TC  4 

104  8(3,3) =-B( 2,2) 

SR(3,1)=-C(5!*B(1,3)  »  ISR(3,1)=0 
SR  (  3 , 1 )  =R  (  2  »  3  I  i  I  SR  (  2  ,  1  )  *0 

SK( 1,1 )=R( 3,3)  %  I  SR ( 1 , 1 ) *0 
CALL  SC  AL  F  ISR(l,l),FSR(l,in 
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call  scale  i sr ( 2 » i ) » i sr  <  2 , i ) ) 

CALL  SCALE  (SR(3,1),ISR<3,1)) 

GO  TO  11 

6  $R(NZ,1)=C(N)-H(  1,3)  $  lSR(N2,i)=0 

SR ( M2- 1 , 1 ) =-C ( N1 > -B { 2  » 3  I  $  I  SR  f N2-1 , 1 ) *0 
CALL  SCALE  { SR ( N2 , 1 ) , I  SR ( N2 , 1 )  ) 

CALL  SCALE  (SR(N2-1,1),ISR(N2-1,1>) 

IF  (N2-2)  11, 11, 9 

9  00  10  J=3,N2 
K*N2*1-J 

SRIK.l >«-R(J,3)  »  I  SR ( K  » 1 ) *0 
CALL  SC  ALE  < SR ( K , 1 ) , I  SR ( K , 1 ) ) 

10  CONTINUE 

11  RETURN 

12  SR < 1 » 1 ) =C ( 1 )  i  ISRI1(1)*0 
SR ( 2 » 1 ) *-C ( 2  )  t  I  SR ( 2 , 1 ) *0 

13  CALL  SCALE  ( SR ( 1 , 1) , I  SR ( 1 , 1  )  ) 

CALL  SCALE  (SR (2,1), I  SR (2,1)) 

GO  TO  11 

14  SR ( 1 , 1 ) *-C ( 4 ) 

SR  (  2, 1 ) *C ( 3 )-C  ( 1 ) 

I  SR ( 1 , 1 ) *0 
t  SR ( 2 , 1 ) *0 
GC  TC  13 
END 

SUBROUTINE  RECON  ( ROOTS ,A,IA,0*N) 
DIMENSION  R00TS(2.50) ,D( 51) 

X*A1I X* I  A 

CALL  UNSCAIE  (X.IX) 

DO  i  I  * i  ,  N 
0 (  I)=C.C 

1  CONTINUE 
UlN  +  l ) *1 «  0 
I«1 

NL»N-i 

2  IF  ( ROO I S ( 2 , I ) )  3,7,3 

3  r«ROCTS(l,n*ROOTS(l,I) 

U*KOO  T  S ( ? , I ) *R00  TS ( 2 , I ) 

T«T»U 

U«2.0*ROOTS( i  ,  l  ) 

UO  5  J* i *  NL 

IF  ( I ♦ J-N )  5,4,4 

4  U( J  » -0 ( J*2)*T*0( J) 

O' J)«0( J)-U*0( J»l I 

5  CONTINUE 

0( N ) *  T*0 ( N ) 

0(N)«0(N)-U*0(N»i) 

0(N*1 )»T*0(N*1 » 
l«l  *2 

6  IF  ( N-  I  )  10,2,2 

7  00  9  J* i  ,  N 

IF  ( J* I -N )  9, 9, A 
3  0(J)*0(J*l)-D(J)*ROOTS(l,I) 

9  CONTINUE 

OIN*l)«-OfN*l)*ROOTS(t»l » 

I  ■  I  »4 
GC  ro  6 

10  NS«N»l 

00  11  I  I « 1 , NS 
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D<  t I )=D(  II  )*X 

11  C  C  N  T I N  U  C 
«F  FUKN 
END 

SUPRCUT  I  NF  UUADIV  (  N,  A  ,  I  A  ,  R  , 0 , 0  ) 

D I  MENS  I  IN  A(51,?>,lA(5i,3),R(2),D(2),ft<4)) 

0  (  1  )  =  A  (  l.miRM  A(1  ,1  ) 

CAI.L  UNSCALE  (R(l), IB) 

IF  ( M  —  2  )  4,4,1 

1  AA=A{2,i)UAA=IA(2,l) 

CALL  UNSCALF  (AA,IAA) 
rt( 2 ) = AA-R ( 1 )*D(  1  I 

IF  (N-3)  4,4,2 

2  NT  =  N- 1 

DC  3  I -ot NT 
XN  =  R I  I -i 1*0(1) 

YN  =  H  (  1-2 )*D(2  » 

AA=  A (  I , L ) i I AA -  I  A ( 1,1) 

CALL  UNSCALE  ( A A , I  A  A ) 

9  (  I  ) =AA-( XN+YN) 

3  CONTINUE 

4  XN  =  fl(N-i  )  *0 (  1  ) 

YN  =  H  ( N-2 ) *D( 2  ) 

A  A  =  A  (  N , 1 ) 1 1  A  A  *  I  A (  N »  1  ) 
call  UNSCALE  (AA.IAA) 

R ( 1  I =AA- ( XN4YN) 

AA  =  A(N*1,  Dll  AA=  I  A(  N*1 ,1  ) 

CALL  UNSCALF  UA.IAA) 

R(2)*AA-MN-1 ) *0  I  2 ) 

RETURN 

END 

SUPRCUTINF  UCURLCG  (X.IX.Y.IY) 

r  *  64  •  0 

IF  (  X  )  1 , 5 , 3 

1  PRINT  2 

2  FORMAT  U6HCTHE  LOG  OF  A  NON-POSITIVE  NUHBFK  IS  REQUESTED) 

5  Y  *  0 . 0 
I  Y  -  0 

GO  TC  4 

3  TC*IXl 

Y»AlDG(  X)»TO*ALOr,(T) 

I  Y*0 

CALL  SCALE  (Y,IY) 

4  RETURN 
FNO 

SUPROUTINF  DCURLFXP  IX, IX, 2,12) 
i  E  XPI X ) *| 2*« 

IF  (IX!  I, A, 6 

1  I *- l Xt I  1*6*1 
DO  4  J*  i  »  I  I 

X  ■  X  M  0  D  F ( 12,2) 

IF  (X)  2,3,2 

2  12*12-1 
2*64.0*2 

3  12*12/2 
2*DiORTI2) 

CALL  SCALE  (2,12) 
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4  con  r  i  Nut- 

5  RETURN 

6  1  I  X 

00  7  J*l,l 
Z«Z*ZilZMZ*IZ 
CALL  SCALE  IZ,IZ> 

7  CONTINUE 
GO  0  5 

8  CALL  SCALP  (Z,IZ) 

GC  TO  5 

END 

SUBROUTINE  ADD  ( X , I  X, Y, I Y  ,  Z  .  I Z > 

IF  (X)  3,1,3 

1  Z«Y 

I Z  *  I V 

2  RETURN 

3  IF  m  5,4,5 

4  Z«X 
IZ«IX 
GC  TO  2 

5  IDIFf-IX-lY 

IF  UOIFF)  6,7,7 

6  I  A*  I  Y 
A»Y 
8*  X 

IDIFF*-I0IFF 
GO  TC  8 

7  I  A- I  X 
A-X 
B*Y 

8  IF  U6-I0IFF)  9,9,10 

9  Z« A 

I Z  *  I  A 

GC  TO  2 

10  IF! 1 0 1 F  F )  11,12,11 

11  DC  12  tM.IDIFF 
8*8/64.0 

12  CONTINUE 
Z*A*ft 

IZ  ■  I A 

CALL  SCALE  U,!Z> 

GO  TO  2 
ENO 

SUBROUTINE  SUBTRACT  ( X , I  X , Y , I Y , Z ,  I  Z  I 
y-Y 

CALL  ADO  (X, IX, fa, IV, 2,12) 

RE  TURN 
ENO 


subrout i nc  scale  «x,m 

TYRE  DOUBLE  X,Y 
REC64*)  :/64.C 

IF  (X)  1.11,2 

1  Y  •  -  X 

GO  TC  3 

2  Y* X 

3  IF  (64.0-Y)  4,5,5 

4  Y-V/64.C 
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GO  TO  3 

5  (F  ( Y-RCC64)  6,7,7 

6  Y 1 Y*  64 . C 
I  X  =  l X-  1 

GO  TO  5 

7  IF  IX)  3,9.  ) 

3  X  a  -  Y 

GO  n  L 
9  Xs  Y 

10  Rfc  TURN 

11  I  X  3  G 

go  ro  n 

FNQ 

SUf'RCl'TI*Jr  UN6CME  IX,1  X) 

IF  ( I  X ♦ 3 4  )  3,4,4 

3  XaC.O 
I  X  *  3 

GO  TO  ? 

4  IF  ( IX-S4  )  5,5,6 

6  X=1 . OF ♦ 1 53 
I  X*0 

PRl.I  7 

7  FORMAT  (25H3EXP.  OVERFLOW  IN  UNSCM-fc) 
GO  TC  0 

5  IF  l  I  X )  1,2,1 

1  X*X*64. 0 *•  I  X 

3  l  X  =  0 

2  RETURN 
ENO 


APPENDIX  B 


C-907  FINDING  ZEROS  OF  REAL  POLYNOMIALS 
MAIN  PR  oG  RAM 

DOES  INPUT  AND  .'UTPUT,  AND  CALLS  RSSR 

DIMENSION  A( 51,3) , I A( 51 ,3 ).ROOTS( Z. 50) , T ITLfc ( 10) ,DI El ) 

TYPE  INTEGER  DEGREE 

1  READ  7,  TITLE 

2  FORMAT  liOAS) 

4  PRINT  2,  TITLE 

5  FORMAT  UH1,1CAD) 

READ  6  j  DEGREE, .M,MMAX, DELTA, EPSILON 

6  FORMAT  (?I6,2tl2,4$ 

IF  ( M )  7,7,8 

7  M  =  10 

MM A  X  =  1 5 
DEL  TA  =  0  #  OOC  » 

EPSILON =0,000001 

8  PRINT  9,  DEGREE, M.MMAX, GEL TA, EPSILON 

9  FCKMAT  (  AH0DCGREE=I3,17l!  NO.  OF  MOOT  SU=I3,7H  MMAX=*7,/ 
17H  DELTA  =  E12. 3.1CH  EPS  I LON  =  E 1 2. 3  ) 

Nl=OFGREF+l 
PRINT  1C 

10  FORMAT  ( i  9  H  D I  N  P  U  T  COEFFICIENTS) 

12  DO  14  1=1, N1 
READ  12,  X 

12  FORMAT  lEie.IG! 

A  (  I  ,  1  )  =  X 

l A ( I,1)=0 
J  =  I-1 

PRINT  13,  J , X 

13  FORMAT  (3HCAI I2,3H)=  1PE18.10) 

CALL  SCALE  (A(I,1),IA(I,1)) 

14  CONTINUE 

19  N2=DEGREF 

CALL  RSSR  ( A, IA.RCOTS, DEGREE, M,MMAX, DELTA, EPSILON, D) 

20  IF  (Nl- (DEGREE  +  1  )  )  29,29,21 

21  PRINT  22 

22  FORMAT  (24HCR00TS  OF  THE  POLYNOMIAL) 

PRINT  23 

23  FORMAT  I42H0REAL  PART  IMAGINARY  PART) 

K=0EGREE*1 

DO  25  I =K , N2 
X=ROOTS( 1, I ) 

Y=RQGTS! 2, I ) 

-—PRINT  24 ,  X  ,  Y 

24  FORMAT  (1H01PE24.15,1PE26.15) 

25  CONTINUE 

IF  (DEGREF)  26,26,2 9 

26  PRINT  1 7 

27  FORMAT  ( 2 5H0R ECO NST RUC T E D  POLYNOMIAL) 

DO  23  1*1 , N 1 

X  =  0<  I  ) 

J  *  I  - 1 

PRINT  13,  J , X 

28  CONTINUE 
GO  TO  1 

29  PRINT  3C 

30  FORMAT  (21H0S0MF  ROOTS  NOT  FOUND) 

GO  TO  1 

End 

SUBROUTINE  RSSR  ( A, I  A, ROOTS, DEGREE, N, MM AX , DEL  T  A, EPS  I LON, D ) 
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DIMENSION  A<  51, 1),  I  A(  51, 3  »,RCC  TS! 2,501  ,lj!51 1  ,RCMCD(5C  )  , KRGMOb!  5 
1NONRTI50) »MNCNRT 150) 

TYPE  INTEGER  DEGREE 
N=DEGRE  C 
IE  !N)  1,1,3 

1  OEGRFC*NCUR 

2  RETURN 

3  N1=N*1 
N2=N1*1 

00  5  l  =  *. ,  N 
K  *  N  a  -  I 

IF  IA|K,il!  6,4,6 

4  J*N1-I 

ROOTS! i , J) =  C.G 
ROOTS! 2, J)=0,0 

5  CONTINUE 
DEGREES 
GO  TO  2 

6  Ni  =K 
N=K-1 
NCUR*N 
NL*N 

7  CALL  ROCTSO  ! A, l A.NCUR, M) 

CALL  REALROGT  ( A , t A , M ,NCUR» Of L T A, E PS I LCN , KOMCO , MRCMJO, 
INONRT.MNONRT.NCC , ROOTS) 

IF  { NCO )  1 2 » i 2  ,R 
S  N1*NCUR*1 

CALL  CCMPROOT  !  A  ,  I  A , RCMOO  ,RCC TS , M, MNONRT , NO NRT , 
i MR OMOD, NCO, DELTA, EPSILON, NC UR  I 
IF  INCUR)  12,12,9 
9  IF  I NL-NCUR )  11,11,10 

10  NL  =  NCl)R 
GO  TO  7 

11  M»M  <■  1 

IF  IMMAX-M)  1,7,7 

12  CALL  RECON  ( ROC TS, A ( 1 , 1) , I A ( 1 , 1 ) , 0, DEGREE ) 

GO  TO  1 

ENO 

SUBROUTINE  ROCTSO  « A , l A , NCUR, MM ) 

DIMENSION  A(51,3),IA(51,3) 

N1*NCUR*1 
00  1  J  *  1 ,  N 1 
A! J,2)«A{ J,i ) 

IA(J,2)sIA(J,i) 

A(J,3)*0«0 
I A !  J ,  3  )  =  C 

1  CONTINUE 

OC  9  M* 1 , MM 
DC  6  J* 1 , N 1 
K1*N1-J 
K2*  Jl-l 

KM*XMIN0F<Ki,K2) 

I F ! KM )  11,4,11 

11  00  4  L* 1 , KM 

LR*XMOOF I L , 2 ) 

JL* J-L 

v)LP  =  J*L 

IF  <LR)  2,2,3 

2  X*A! JL,2)*A! JLP,2) 


* 
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IX=IA(JL,?)+IA(  J  L  P , 2 ) 

CALL  SCALP  !X,!X) 

CALL  AOu  (  A<  J,  3)  ,  l  A!  J,  3  )  ,  X,  I  X,  4!  J,3  )  ,  l  4(  J,-1  )  ) 

GO  TO  4 

3  X=A  <  -IL * 2  ) *A  C  JLP.2) 

I  X= I  A ( J L  »  2 )  ♦ 1  A ( JLP,  2  ) 

CALL  SCAL?  ( X  ,  ! X  ) 

CALL  SUBTRACT  {  A  (  J  ,  3  )  ,  I  A  (  J  ,  3  )  ,  X  ,  I  X  ,  A  (  J  ,  3  )  ,  I  A  (  J  ,  }. )  > 

4  CONTPHH 

A( J, 3) =2.0*A<  J,  ) 

CALL  SCALF  l  A<  J,  3)  ,  I A  (  J,  3) ) 

X=A ( J , 2 ) **2 
IX*IAU,?)MA(J,?> 

CALL  SCALF  (Xf!X) 

CALL  ADO  (A! J,3)  ,IA(J,3)  ,X, IX, A! J,3>, I  A ( J.3)  ) 

Jk=  Xf’CDF  (  J  .  2  ) 

IF  (Jk)  .-,5,6 

5  41J,3)=-A(J,3» 

6  CONTINUE 

IF  (VM-M)  9,9,7 

7  DC  3  J  =  1 , N 1 
A< J , 2  >  =  A C J,3) 

IA( J,2)-IA{ J, j) 

A  <  j  ,  3  )  = : .  0 

I A I J, 3 ) =0 

8  CONTINUE 

9  CONTINUE 
RE  TURN 
END 

SUBROUTINE  RE ALROOT  ( A, I  A , M , NCUR , UELT A , E PS i cCN, ROMOD, MR OMGu 
INONRT.MNDNRT.NCO, ROOTS) 

DIMENSION  A ( 51,3 ) , I A ( 51 ,3) , ROOTS  I  2, 30 ) , ROMOD I  50 ) , MRCMGO! 5» ) 
1NCNRTI50) fMNONRTl 50), RAT  10(51) , l PI  VI 51 ) , ARFu! 50 ) , I ARFOl >0) 
RATIO! I ) =1.0 
DC  4  1=2, NCUR 

I I  =  XMCOF 11,2) 

IF  ( A( I  ,3)  )  2,1,2 

1  RATIO! I  )  =0 . 0 
GO  TO  4 

2  T  =  A ( 1 , 2 ) * A ( 1,2) 

IT*I  AM  ,2  )  ♦!  Al  1 , 2) 

CALL  SCALE  (T,IT) 

T  —  T / A ( 1,3) 

I T»I T- I  A ( ! ,3) 

IF  (IT-2)  50,50,1 

50  IF  ( I  T*2 )  1.51,51 

51  CALL  UNSCALE  ( T, IT) 

RATIO! I ) =  T 

IE  Ml)  3,3,4 

3  RATIO! I ) =-RAT 1 0 1 C  ) 

4  CONTINUE 

RAT  ICINCURM  )  =  1. 0 
I P i V ! 1 ) =1 
I P l V ( NCUR* 1  )  =  1 
DO  7  1*2. NCUR 
X=ABSF (RATIO! I ) - 1 . 0 ) 

IF  (X-DELTA)  5,6,6 

5  IPIVM  )  =1 
GO  TO  7 
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6  ipivh)=c 

7  CON  * ’ NUt 
NCUR1  JCUR+1 
11*0 
MUL 1=0 
1*1 
14*1 

8  I1=I1*1  \ 

1 2* 1 1+ l  1 

MUL  T  =  MUL T*  1 

IF  { I P I V  (  1 2  ) 1  3,8,9 

9  ROMODf |4)=A(I2,3)/A(I,3) 

IR0M0n*lA{  l  2, 3 )- I A(  1,3) 

CALL  SCALF  (RCM00M4)  ,  IROMOD) 

IF  (RCMCDU4))  10,11,11 

10  RGMGIJ (  U)=-RCM0U(I4) 

11  CALL  DOURLOG  IRC MOD ( I  4) , I ROMOO ,XN, I XN) 

T*2  **M* 1 1 

XN=XN/T 

CALL  SCALE  <XN,IXN) 

CALL  OCURLEXP  I XN, I XN , ROMOO ( I  4 ) , S RCMOD ) 

I F { IROMCO-74  )  ICO, IOC, 101 

100  I F ( I RCMgO+74  )  101,101,102 

101  ROMOO 1  I  4 ) *0. C 
IROMOD*C 

GC  TO  1C3 

102  CALL  UNSCALE  ( ROMOO ( I  4 ) , I ROMOO ) 

103  MRCMOD (14) =MUL  T 

IF  (NCUR+1-12)  13,13,12 

12  1=12 

1 4*  I  4-f  1 
MUL  T *C 
11*0 
GO  TO  8 

13  0*0.0 
NCC*C 

DO  22  1=1,14 
KL*  I  *1-  I 
W*-ROMOD{KL) 

I5=MKOMOO«KL) 

DO  20  J *  1 ,  I  5 
J*J 

14  CALL  TEST  (  A , 1  A , W,Q,NCUR , ROMOO ( KL) , C PS  I  LON, K ) 

IF  (K)  17,17,15 

15  ROOTS! 1 , NCUR) *-W 
R00TS(2,NCUR 1=0.0 
AREU ( 1 ) * A ( 1  , 1 ) 

I ARED! 1 ) * l A ( 1 , 1 ) 

00  16  L*Z ,NCUR 
Y*ARED( L-i ) *W 
I Y« I AR  E  0 ( L  ~ 1 ) 

CALL  SCALE  (Y.IY) 

CALL  SUBTRACT  ( A ( L, 1 ) , I  A { L, 1 ) , Y, I Y, ARfcU! L) , I AR EDI L II 
A ( L , 1 ) * AREO ( L ) 

I A ( L , 1 ) *  I  ARE  (L) 

16  CONTINUE 
GO  TO  19 

17  IF  ! W )  18,21,21 

18  W*-W 

GO  TO  14 
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19  NCUR  =  NC  UP  - 1 
2  0  CONTI NUF 
GO  TO  22 

21  NCC=NC0+1 
N0NRT(NC0)=KL 
MNCNKT ( NCC )  =  I  5*1- J 

22  CONTINUE 
RETURN 
END 

SUBROUTINE  COMPROOT  {  A  ,  I  A  ,ROMO[;,  ROC  TS ,  M  ,  V.NONR  T  ,  NCNR  T  , 
lMROMCU,NCO,DELTA,EPSILO.N,NCUK) 

DIMENSION  A(51,3),IAI51,3),R0M0D(5C)fR0}TS(<:,5C)  ,  $R I  5 1 , 3  I , i SR ( 
1  ),  SRONCUI  EG)  ,  SRC0TSI2 ,50  )  »MNONR  T  (  50  ) ,  NO  N  U(  5 C )  ,  KSROMC  0(  ^C  )  , 
2NSCNRT(49' ,MSNCRT(49) ,MRCMOO( 5C  > ,U<2) ,K( i ) ,H(*9) 

DC  23  1=1, NCC 
JA*NONR  T(  l  ) 

1 1  =  MNLNR  T (  I  ) 

11=11/2 
IF  (ID  1,1,2 
1  11  =  1 

2  IF  (RCMOO(JA)J  3,23,3 

3  Q  =  RCVCD«JM 
DC  22  J=1,I1 

300  CALL  SUBRFS  « A ,  I  A , NCUR , SR , I  SR , U ) 

IF  (NCUR-4)  4, ICC, 5 

100  N  SC  UR  =2 
GO  TO  IC1 

4  N  SC  UR  =  1 
J2  =  l 

GO  TC  6 

5  NSClJR=NCUR-3 

101  J2=NSCUR 

6  LL  =NSCUR ♦  1 

IF  (NSCUR-D  7,7,1 

7  IF  (SRI  1,1 D  6,10,3 

8  X  =  SR( 1,1)  t  I X=I SR( 1,  1  ) 

Y  =  SR ( 2 , 1 )  $  I Y= I  SR ( 2 , 1 ) 

CALL  UNSCALE  ( X, I  X) 

CALL  UNSCALF  ( Y , I Y ) 

SRCOTSI 1,1 »=-Y/X 
NSCUR=C 

GO  TC  11 

9  CALL  ROGTSQ  (  SK , I  SR , NSC  UR  ,M  ) 

CALL  RE ALR°  r  (  SR ,  I  SR  ,M  , NSCUR ,  )FLT A , EPS  I  LON ,  SRC  MOD,  V SC f,C  J , 
1NSCNR  T , RSNOR  < , NSCC , SRCO T S I 
IF  (J2-NSCUK)  10,10,11 

10  SROOTSI i,J2)=C.O 

11  SROOTSI  l,J2)  =  SROOTS(l,J2)*RjMOD( JA) 

T  =RCMCD  (  J  A  )  ♦  KOMOC'I  JA) 

IF  (SROOTSI i,J2)-T)  12,16,16 

12  W*SK00TS 1 1 »  J2  ) 

WE*RCMCU(JA)*ROPOO(JAI 

CALL  TEST  » A , I  A , W, WF , NCUR , ROMO D( J A ) , EPS  I L ON , K  ) 

IF  (K>  16,15,13 

13  ROCTSI 1,NCUK)»>W/2.C 
T=4.0*WF 

U*W*H 
T*  T-U 

IF  (T)  201,20 1,202 


* 


i. 
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201  T*-T 
U’DSORT  (T) 

R0CTS(1  ♦NCUK)=R^C  TS(1  ,  NCIJR ) -IJ/ 2  .  G 
ROOTS (l,NCUR-i)=-(W-U)/2.0 
ROC  TS I  2 , NCUk ) =0 . 0 
ROOTS  ( 2  ,  NCI  JR- 1 )  =0. 0 
GO  TC  20* 

202  UxDSQRT  <  T ) 

203  R0CTS(2,NCUk)*U/2.0 

ROOTS ( 1  *  NCUR- I )  =R00TS 1 1 » NCUR) 
ROOTS(2,NCIlR-i)*-ROCTSC2,NCUR) 

204  0(1)*W 
0«2)=WE 

CALL  QUAf)  f  V  INCUR, A,  IA,R,0,R) 

JX*NCUR-1 
00  14  JY*1 • JX 
A( JY»1 ) *B( JY) 

I A  (  J  Y .  1 »  -  C 

CALL  SCALE  ( A ( JY , 1 ) , I  A ( J Y  ,  1  ) ) 

14  CONTINUE 
NCUR*NCUR-2 
GO  TO  22 

15  W*-W 

CALL  TEST  ( A , I  A , W , WE , NCUR .RGMOUl  JA ) , EPS  I  LON, K ) 

IF  (K)  16,16*13 

16  IF  { J2-INSCURU)  )  23,17,10 

17  IF  (J?-l)  23,23,13 

18  J2* J2-1 
SKCOTSI 1 • J2 ) *0,0 
GO  TO  12 

10  IF  (SROCTSUf  J2)-SR00TS(1,J2-1  M  2C,21,20 

20  J2* J2-1 
GO  TO  11 

21  J2* J2-1 
GO  TO  16 

22  CONTINUE 

23  CONTINUE 
RETURN 
END 

SURRCUTINE  TEST  ( A, I  A ,W , Q,N,R3 MOO, E PS  1 LON,K ) 

DIMENSION  A151,3),IAI51,3),0(3I,IM(3I,TI2I,EI2I,C(51) 
BI1)*C»C4IX«0HW«0 
1 8 ( 11*0 
R<2)*AIi,l> 

IR(2)*IA(1,1) 

00  2  I  « 1 ,  N 
X*W*B  f  2 ) 

I X* I R ( 2 ) 

CALL  SCALE  <X,|X) 

Y»0*P( 1 ) 

I Y» I B I  1 1 

CALL  SCALE  IY,IY) 

CALL  ADD  (X,IX,Y,|Y,Z,IZ) 

CALL  SUBTRACT  I  A ( I *1 , 1 ) , I  All ♦ 1 , 1) , i , 1 1 , B I  3 ) ,  I C. < 3 >) 

IF  (N-l)  2,2,. 

1  8<il*R<2l 
I B  (  1 )  *  1 0  I  2  ) 

BI2  )*6< 3) 

IBI2I-IBI3I 
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2  CONTINUE 

3  KCUNT=1 

CEPSl l=CPSl LGN 
T  ( 1 ) =C«  Gt T ( t ) =0. C 
N1  =  N  +  1 

X=2.0*EPS!LGM 

Y=X*RCMGD 

E  (  1 )  =RCM3D*Y 

E  (  2  )  =R0M0r»CEPSIL*R0M0D 

0G  7  1=1, Ml 

IF  (  A  (  I  ,  1  )  )  4,5,5 

4  Cm=-A(I,i)ilC=IA(l,n 
GC  TO  6 

5  C(  l ) =A( I , 1 ) 1 1 C = I  A (I » 1  ) 

6  CALL  UNSCALF  (Cl  I) , IC  ) 

T(i)*T(i)*E(l)+C(l) 

T( 2 ) =T( 2 >*E( 2  )  +C  (  I  ) 

7  CONTINUE 
0IF=T(1)-T<2> 

3  If  (G)  15, 9, IE 
9  IF  (  P  (  3  ) )  10,11,11 

10  B  (  3  ) =-B I  3  ) 

11  IF  ( I R ( 3 ) ~74  )  ICO, IOC, 12 

100  IF  ( I P ( 31+74  I  12,101,101 

101  CALL  UNSCALE  (8(31,16(3)) 

IF  (DIF-RUM  12,12,17 

12  K  =  0 

IF  (K0UNT-2)  *3,16,16 

13  IF  (0)  **,14,15 

14  IF  (W)  l  5 , 1 6 ,  j.  6 

15  SENSE  LIGHT  2 
K0UNT=KCUNT*1 

16  RETURN 

17  K  =  1 

GO  TC  lo 

18  IF  IIPI2I-74  )  102,102, *2 

102  IF  ( I B ( 2 ) *74  )  12,103,103 

103  CALL  UNSCALE  (8(2),IB(2)> 

IF  I l R I  31-74  )  104,104,12 

104  IF  (IB(3J*74  )  12, 105, 115 

105  CALL  UNSCALF  18(31,10(3)1 
X«0*0(2)*R(2> 

Y  =  W  *  P  <  2 )  *  B  ( 3) 

2  =H ( 3  )  *0 ( 3  ) 

Y*X-Y*2 

IF  (Y)  19,17,20 

19  Y*-Y 

20  D i F  *DI F  *Ol F 

IF  (OiF-Y)  12,17,17 
ENO 

SUBRCUT i Nc  SUBRFS  (  A  ,  I  A , N , SR , I  SR , KC MO 0  ) 

0 1  MENS  I  jN  A  (51 , 3  ) ,  I  A(  51  ,  3)  ,  SR(  51 ,3  )  ,  I  SR  (  51, 3  )  ,  C  (  5  l  )  ,B(  5-.  ,  3  ) 
Ni 1 
T  *  l  .  C 

00  l  I  *  » , N 
J*N 1- i 
T  *  T  *RGMuP 
C( J)=A< J,1)*T 
IC*IA( J,l) 
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CALL  UNsCALE  (C(J),IC) 

1  CONTINUE 
C(Nl)*A(M,l) 

IC  *  I A  (  N 1 « 1 ) 

CALL  UNSCALE  <C(Ni),lC) 

IF  IN-2)  12,  12,  2 

2  N2*N-2 

DO  3  |*1, N2 
81 l , 1 ) *C.  3 
8(1,21*0.0 

3  CONTINUE 
1*2 

B(l,2)*Cm 
A  B( 1 , 3 ) *C | I ) -B 1 1 » 1 ) 

00  S  J*2,N2 

8(J,3l*-8(J-l,2l-B(J,l| 

5  CONTINUE 

IE  (N-(3+m  IOC, 6, 6 

6  1*1+1 

DO  7  J*1,N2 
8(J,1J*R(j,2) 

8(  J,2)*BU,3) 

7  CONTINUE 
GO  TC  4 

100  IF  IN- 4)  14, 101, Q 

101  IF  ( N- ( 2 ♦ I ) )  104,102,102 

102  1*1+1 

OC  103  J*i ,2 
8(J,1)*H(J,2) 

HI J,2)*B( J, 3) 

103  CON T I NUt 
GO  TO  4 

104  6(3,3) *-B ( 2,2) 

SR(3»l)*-C(5)+8( 1,3)  X  IS«I3,1)*3 
SRI  2,1 )*B(2,3 )  X  ISK(2,1)«C 

SRI 1«1)*B(3,3)  \  I  SR  I  1 , 1 ) *Q 
CALL  SCALE  I SR( 1,1), I  SRI  1,1)) 

CALL  SCALE  I  SR  I  2 , l ) , I s R I  2, 1 ) ) 

CALL  SCALF  I  SRI  3,1 1 , ISK( 3, 1 ) ) 

GO  TO  li 

8  SR ( N2 , 1 ) «C I N) -fl ( I , 3 )  X  (SR(N2,1)*C 
SRIN2-l,i)**C(Nl 1-8(2, 3)  »  ISR(N2-i 
CALL  SCALE  (SR(N2,1 I, ISR(N2,1) ) 

CALL  SCALE  ( SR(N2-1,1 ), l SHIN2-1 ,1 ) ) 
IE  (N2-2)  11,11,9 

9  DC  1C  J»3,N2 
K-N2  1-J 

SRU.l)— R(  J,  J)  $  I  SR  (K  ,  1  )«0 
CALL  SCALE  ISR(k,1),|SR(k,1)) 

10  CCNTINUC 

11  RETURN 

12  SKI  1,1 ) *CI 1 )  %  I SR< l, 1 ) *0 
SKI2,1 )«-CI2)  X  I  SR ( 2 , 1 ) *0 

13  CALL  SCALE  I  SR  1 1 , l ) , I  SR ( 1 , 1  ) ) 

CALL  SCALF  I  SR  I ?  ,  1 ) , I  SR  I  2 , 1 ) ) 

GO  TO  11 

1*  SR(»,1)*-C(4) 

SRI  2, 1 ) *C I  3 ) -C 1 1 ) 

I  SR  1 1 , 1 ) *0 

I  SRI 2,1 1 *0 


,  1 )  *  j 
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GC  TC  13 
END 


SURRCUT I NF  KFCON  ( ROC TS , A , l A, o , N ) 

UIMLNSION  RCGTSC2.50) ,0(51) 

X  =  A  4  I  X  =  I  A 

CALL  UNSCALL  ( X, IX) 

OC  1  1*1, N 

0 (  U=c.w 

1  CONTINUE 

o  (  n  + 1 )=i, : 

1  =  1 

NL  =N- i 

2  IF  CP00TSI2.I I)  7 , 7 , 3 

3  T«RCCTSC:,I)*RCCTSC,II 
U=RCOTS( I  )  *R00  T  S ( 2  , I ) 
t  =  t+u 

U  =  2.0*RCCTSC  ,1  ) 

OC  5  J  *  «.  »  NL 
IF  (I+J-N)  5,^,4 

4  OC  J) =D( J+? ) +  T*D( J) 

0( J)=D( J)-U*D( J+l) 

5  CONTINUE 

0  (  N  )  =  T  *  0  (  N  ) 
n<N)=0(N)-U*U(N+l ) 

0  (  N  ♦  1 ) »  T  *D ( N* 1  ) 

1  =  1+2 

6  IF  (N-I)  10,2,0 

7  00  9  J  *  1  ,  N 

IF  ( J ♦ I  - N )  9,N,“ 

8  0( J) =C( J+l ) -0 ( J) *K0CTS( 1 ,  I  ) 

9  CONTINUr 

0 (  N ♦  1  )=-l>(N+l  )  *R0OTS<  1,1) 

1  =  1+1 
GC  T  C  6 

10  N  S  =  N  ♦  1 

OC  11  11=1, NS 
IH  I  I ) =0 ( 1  I ) • X 

11  continue 
RETURN 
cNC 

SOPKOUTINF  w  1 1 A  0  I  V  (N, A, I A,R,0,P) 

01  HENS  I  ON  A(5i,3),|A(5l,3),R(2),u(2),P(49) 
8(1  )  =  A  *  1,1)418*1  AU, 1  ) 

CALL  UNSCALE  18(1), IH) 

IF  (N-2)  4,4,1 

1  AA-AC  2, 4 ) i I  A A ■  I A  I  2 , 1 ) 

CALL  UNSCALf  (AA.IAA) 

9(2  )  *AA-«M  1  )•()(!  ) 

IF  ( N- 3  1  4,4,2 

2  NT*N- 1 

OC  3  I  *3, NT 
XN*  8 (  1  —  1  )  *0 (11 
YN«B( I  - 2 ) *0 ( 2  I 
A  A  *  A (  1,1  )t|AA  =  ]A(  1,1) 

CALL  UNSCALf  ( AA , I  A A ) 

R(  I  ) *  AA- ( XN  +  YN) 

3  CONTI  NUf 

4  XN*8(N- 1 )*D( 1 > 


* 


i 


ii 


YN*RlN-2 )*0(?) 

AA*  A  (  M  *  1 ) if  AA* I  A ( M» 1 ) 

CALL  UNSCALE  ( A A ( I  A A ) 

R(1)*AA-IXN*YNJ 

AA*  A  (  N*  1 , 1 ) S I AA  = | A ( N ♦ 1  *  1  ) 

CALL  UNSCALE  (  AA  ,  I  A  A  ) 

R(2)»AA-R(N-1 )*D(2> 

RETURN 

ENo 

SUPR9UT l NE  OCUBLCG  CX,IX,Y,[Y) 

T*64.C 

IF  I  X  )  l ,  5 , 3 

1  PRINT  2 

2  FORMAT  (4BHCTht  LCG  OF  A  NCN-PC S I T 1 VF  NUMOFR  IS  RFQIIE  j  Thi> ) 
5  Y*C . 2 

IY«C 
GC  TO  4 

TO* I Xi 

Y*ALCG(X)*T04ALOG(T) 

I  Y* C 

CALL  SCALE  *  Y ,  I  Y  ) 

4  RETURN 
ENU 

SUBROUTINE  UCUBLEXP  (X,I  X,Z,IZ) 

z*  Expixmz*c 

IF  (IX)  l.P.6 

1  1 *- I Xi I  1*6*1 
DC  4  J*  I ,  I  I 
K«  XHODF (  I Z  «  2  ) 

IF  (K)  2,3,2 

2  IZ-IZ-l 

l *64, C*Z 

3  IZ-lZ/2 
Z*OSCRT ( Z I 

CALL  SCALE  <z,m 

4  CONTINUE 

5  RETURN 

6  t»fc*IX 

oc  7  j«;,i 
2*f*ZftlZ*l2«IZ 
call  SCALE  (Z.IZ) 

7  CONTINUE 
GO  TC  5 

«  call  SCALE  (Z,IZ) 

GC  TC  5 
ENO 

SUBROUTINE  AUO  I  X , I  X , Y , I Y ,Z , 1 Z ) 

IF  (X)  3,1, 3 

1  Z«Y 

I  Z  *  1  Y 

2  RETURN 

3  IF  (Y)  5,4, S 

4  Z»X 
IZ-IX 
GG  TC  2 

5  I05FF-IX-IY 

IF  (IOIFF)  f,7,7 
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6  I  A  -  I  Y 
A  =  Y 

n~x 

i  n  i  f  f  -  i  p  r  f  ► 

gg  rc  p 

7  I  A 1 1  x 
A  =  X 

ft  =  Y 

3  U  (  lt>—  ItJlFFl  9,9,  K 
9  Z  =  A 

I  /  =  I  A 
GC  TO  2 

10  I F ( I D I F  F  I  11,12,11 

11  OC  12  I  -  1  »  I  U  I  E  F 
R=P/64. 0 

12  CONTI NUl 
l- A  +  « 

I I  =  I  A 

CALI  SC  alt  I  2  ,  I  Z  ) 

gc  rc  : 

END 

SUP  R  CU  T  INF  SUBTRACT  (  X  ,  l  X  ,  Y  ,  I  Y  ,  l  ,  l  Z  I 
*  =  -  Y 

CALL  ado  (  X,  t  X  ,W,  |  Y  ,1  ,1 l  ) 

RETURN 

ENC 

SUBROUTINE  AflL[  (X.IX) 

TYPr  CuU«LE  X ,  Y 
RE C  t4  =  1  .9/6-.0 
IE  (XI  1,11,2 

1  Y  =  -  X 

GC  TO  1 

2  Y-  X 

3  IE  (oF.C-Y)  4,5,5 

4  Y  =  Y  /  6  4  .  . 

I  *  =  I  <  ♦  1 

GC  I  '  x 

5  IE  IY-R/C64J  6,7,7 

6  Y  *  Y  •  A>4  ,  . 

1  X  *  I  X  -  1 

GC  TC  5 

7  IE  (  X  )  ft  ,  '•  ,  7 
a  x  *  -  v 

gc  rc  lc 

9  X  *  Y 

10  Rf TURN 

11  I  X  =  0 

GO  T  ^  1  ^ 

enu 


SUP KOUT  IN.-  UNSCAIF  (X,IX) 
IF  (  I  X ♦ H  4  |  3,4,4 

3  X*C.O 

I  X  =  v 

■1C  Tj  2 

4  IF  ( IX-64  »  5,5,6 

6  X» 1 ,CE ♦ l 53 

IxO 
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PRINT  7 

7  FORMAT  (251-UtXP.  OVERFLOW  IN  UNSCALE) 
GO  TO  2 

5  IF  (IX)  1,2,1 

1  X»X*6*.0**I  X 

8  IX  =  0 

2  RETURN 
END 


appendix 


INPUT 

EXAMPLE  1 

6  2C  2C  l.OOOE-4  1. COOE-6 
7056. 

-51072. 

22559.36 

661360.64 

-2101084.2 

2490368. 

-1048576. 

OUTPUT 

EXAMPLE  1 

DEGREE2  6  NO.  OF  ROOT  SQ=  20  MMAX*  20 
DELTA2  0.100E-03  EPSILON2  0.100E-05 

INPUT  COEFFICIENTS 

A (  0)=  7.C56000CCCOF  03 

A (  1)=  -5. 1072COCCCCE  04 

A (  2  )  -  2.255936CCCCE  04 

A(  3  )  =  6.6136C64CCCE  05 

A (  4 ) =  -2. 1C1C842GGCE  06 

A (  5 ) =  2 . 49C 36 0COCOF  06 

A (  6 ) 3  -1.048576000CF  06 

ROOTS  OF  THE  POLYNOMIAL 
REAL  PART 

-3. 99968751 C094593E  00 
3 . 99874995 73142Q6E  00 
2.285716169209813E  00 
2.285716169209813E  00 
1.357365320213065E  00 
1. 3102 36461 141304E  00 


IMAGINARY  PART 

O.OOCOOOOOOOOOOOOE  00 
O.OOOOOOCCOOOOOOOE  00 
O.OOOOOOOOOOCOOOOE  00 
O.OOCOOOOOOOCOOOOE  00 
O.OOOOOOCOOCOOOOOE  00 

O.OOOOOOOOOOCOOOOE  00 
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OUTPUT 


RECONSTRUCTED  POLYNOMIAL 
A(  0)*  7.05609G0CC0E  03 

At  1>*  -5.1072C09377E  04 

A(  2  )  *  2.2559385004?  04 

At  3 1  *  6.6136077332E  05 


At  4 ) »  -2.101CP46001E  06 
At  5)*  2.4903682667E  06 
At  6)*  -1.048576CCG0E  06 


Unclassified 
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It).  ABSTRACT 

Thin  report  describes  the  FORTRAN  Subroutine  POLYR  and  a  related  complete 
program  BRL-RSSR  for  finding  all  roots  (real  and  complex)  of  a  real  pjlynomial 
equation  N 

P(x)  ~  )  A.xN_l  =  0. 

.  i 

i=0 

The  method  which  is  used  combines  the  root  squaring  and  the  subresultant 
(extracting  quadratic  factors)  procedures  for  calculating  all  roots,  even  multiple 
roots,  of  real  polynomials.  Reconstruction  of  the  coefficients  from  the  product 

N 

II  (x  -  Root.)  , 
i=  1  1 

which  is  included,  serves  as  a  means  of  checking. 


The  described  subroutine  and  program  were  adapted  from  a  similar  program  in 
double  precision  described  in  Report  ANL  6987,  Argonne  National  Laboratory 
by  Erwin  H.  Bareiss  and  Ronald  Hamelink. 
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